home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
listings
/
v_09_05
/
9n05085a
< prev
next >
Wrap
Text File
|
1991-03-20
|
12KB
|
269 lines
/**********************************************************
*
* Copyright (c) 1991, John Forkosh. All rights reserved.
* Released to the public domain only for use that is both
* (a) by an individual, and (b) not for profit.
* --------------------------------------------------------
*
* Function: lint1d ( x0,dx,n, f,y )
*
* Purpose: Creates a table containing n values,
* starting with y[0]=f(x0), y[1]=f(x0+dx),
* through y[n-1]=f(x0+(n-1)*dx).
* However, the value in y[i] won't
* necessarily exactly equal f(x0+i*dx)
* as indicated above. Instead, the table
* is constucted to minimize the mean square
* error due to sampling random points in
* f(x)'s domain.
*
* Arguments:
* x0 (I) Double containing the first point
* in the function's domain.
* dx (I) Double containing x-interval
* between points in the table.
* n (I) Int containing the number of
* points in the table.
* f (I) Address of function returning
* double whose values are to be
* represented in the table.
* y (O) Address of double returning the
* table of n values, as discussed.
*
* Returns: (int) 0 for a normal solution, or
* 1 for any error (e.g., a singular
* set of equations).
*
* Notes: o See the accompanying article for a complete
* description of lint1d().
* o Set #defined value of TESTDRIVE to 1/TRUE
* to compile a sample main() (see below).
*
* Source: LINT1D.C
*
* --------------------------------------------------------
* Revision History:
*
* 01/07/91 J.Forkosh Installation
*
*********************************************************/
/* --- standard headers --- */
#include <stdio.h> /* need NULL ptr value */
#include <stdlib.h> /*for malloc() prototype*/
#define msglevel 1 /* to printf debug info */
/* --- set testdrive to compile (or not) test main() --- */
#define TESTDRIVE 1 /* 1=compile it (0=not) */
/* --------------------------------------------------------
Entry point
-------------------------------------------------------- */
int lint1d ( x0,dx,n, f,y ) /* returns 0=ok, 1=error */
double x0,dx; /*1st point and interval*/
int n; /* number of points */
double (*f)(); /*func to be interpolated*/
double *y; /*table for interpolation*/
{
/* --- local allocations and declarations --- */
int iserror = 0; /* no error detected yet */
int i,j; /* row,col indexes */
double x; /* arg for f(x) */
/* --- need temporary array for coefficient matrix --- */
double *a = (double *)malloc((n*n+1)*sizeof(double));
/* --------------------------------------------------------
Set up coefficient matrix as per discussion in article
(since the matrix is symmetric, the columnwise requirement
for simq() input is irrelevant).
-------------------------------------------------------- */
/* --- first check that malloc() was successful --- */
if ( a == NULL ) /*failed to malloc matrix*/
return ( iserror=1 ); /*return error to caller*/
for (i=0;i<n*n;i++) a[i]=0.0; /* init array to zeros */
/* --- 4's on diag (2's at extremes) and 1's offdiag --- */
for ( i=1; i<n; i++ ) /* skip 0,0 element */
{
j = n*i + i; /* index of diag element */
a[j] = 4.0; /* set diagonal element */
a[j-1] = a[j+1] = 1.0; /* and off-diag elements */
} /* --- end-of-for(i=1...n-1) --- */
a[0] = a[n*n-1] = 2.0; /* set 1st,last diag */
a[1] = a[n*n-2] = 1.0; /*and remaining off-diags*/
/* --------------------------------------------------------
Set up vector of constants as per discussion in article
-------------------------------------------------------- */
/* --- interior points --- */
for ( i=1; i<n; i++ ) /* loop skips 1st point */
{
x = x0 + dx*i; /* next arg to be tabled */
y[i] = 2.0*(f(x-0.5*dx)+ /*const vector calculated*/
f(x)+f(x+0.5*dx)); /* as per article */
} /* --- end-of-for(i=1...n-1) --- */
/* --- boundary points --- */
y[0] = f(x0) + 2.0*f(x0+0.5*dx); /* first x in domain */
y[n-1] = f(x) + 2.0*f(x -0.5*dx); /*use last x from loop*/
/* --------------------------------------------------------
Display input to simq (for debugging purposes if necessary)
-------------------------------------------------------- */
if ( msglevel >= 9 /* want debugging output */
&& n < 10 ) /* have room to fit it */
{
printf("lint1d> input to simq...\n"); /* stub info */
for ( i=0; i<n; i++ )
{ /* --- display row (really col) of the matrix --- */
for ( j=0; j<n; j++ ) printf("%6.2lf",a[n*i+j]);
/* --- followed by literal for y* and const y --- */
printf(" y*[%d] %8.2lf\n", i,y[i]); }
} /* --- end-of-if(msglevel>=9) --- */
/* --------------------------------------------------------
Solve the simultaneous equations
-------------------------------------------------------- */
iserror = simq(a,y,n); /* y[] returns solution */
free ( (void *)a ); /*don't need coeff matrix*/
return ( iserror ); /*back to caller with y[]*/
} /* --- end-of-function lint1d() --- */
#if TESTDRIVE
/**********************************************************
*
* Purpose: Test driver for lint1d().
* Prepares a table of values for
* interpolation of the function f(x)=x*x,
* and compares the resulting accuracy
* with usual linear interpolation.
*
* Command-Line Arguments:
* x0 (I) Double containing the first point
* in the function's domain
* (default=-10.0).
* dx (I) Double containing x-interval
* between points in the table
* (default=1.0).
* n (I) Int containing the number of
* points in the table
* (default=21).
* expon (I) Int containing the exponent for
* the test function f(x)=x**expon
* (default=2).
*
*********************************************************/
/* --- program defaults (x0,dx,n,expon from cmdline) --- */
#define X0 (-10.0) /* 1st point in table */
#define DX 1.0 /* table interval */
#define N 21 /*number pts in table*/
static int expon=2; /* test f(x)=x**expon */
#define NMAX 99 /* easier than malloc */
#define NRMS 101 /*pts per dx for rms calc*/
/* --------------------------------------------------------
Entry point
-------------------------------------------------------- */
main ( argc, argv )
int argc;
char *argv[];
{
/* --- local allocations and declarations --- */
double x0=X0, dx=DX; int n=N; /* defaults */
int i,j; /* loop indexes */
double x; /* arg for f() */
double f(); /* test function */
double y[NMAX]; /* table from lint1d() */
int iserror=0; /*return flag from lint1d*/
double frms=0.0,yrms=0.0; /* mean square errors */
double xa,xb; /*interval bounds for rms*/
double sqrarg; /*dummy arg for sqr macro*/
/* --- linear interpolation (u(x=xa)=ya, u(x=xb)=yb) --- */
#define u(x,xa,ya,xb,yb) (ya*(xb-x)+yb*(x-xa))/(xb-xa)
/* --- square a double --- */
#define sqr(x) (sqrarg=(x),sqrarg*sqrarg)
/